Does the continuum theory of dynamic fracture work? 
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We investigate the validity of the Linear Elastic Fracture Mechanics approach to dynamic fracture. 
We first test the predictions in a lattice simulation, using a formula of Eshelby for the time-dependent 
Stress Intensity Factor. Excellent agreement with the theory is found. We then use the same 
method to analyze the experiment of Sharon and Fineberg. The data here is not consistent with 
the theoretical expectation. 
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Recently, there has been renewed interest in the 
physics community concerning the problem of dynamic 
fracture. This interest was kindled by experiments jlj |^| 
showing a universal transition to crack branching and 
the realization that one could not approach this question 
within the confines of traditional fracture mechanics ||. 
Since then, there have been other novel experimental 
findings [Q as well as further evidence of theoretical in- 
adequacies. 

In the standard approach, (see, for example, ||) 
dubbed Linear Elastic Fracture Mechanics, one assumes 
that continuum elasticity is valid everywhere outside of 
an microscopically sized process zone. The instantaneous 
crack tip velocity is then postulated to depend only on 
the singular part of the stress field obtained by solving 
this macroscopic continuum problem. This singularity of 
the stress field is universal in nature, up to an overall 
multiplicative factor, the so-called stress-intensity factor 
(SIF) . Even if there is a well-defined relationship between 
the crack velocity and the SIF, the theory cannot predict 
the form of the relationship, as that depends explicitly 
on physics at the scale of the process zone. As real- 
ized initially by Slepyan||, one way to remedy this defi- 
ciency is to model the entire system as a lattice of mass 
points connected by nonlinear springs. On scales large 
compared to the lattice spacing, the displacement ap- 
proaches that predicted by the continuum theory; on the 
scale of the crack tip, the stress field divergence is reg- 
ularized, thereby allowing for the imposition of a phys- 
ically sensible breaking criterion. This criterion is usu- 
ally in the form of a critical spring displacement, as this 
allows for the possibility of analytical solutions of the 
model@, 1 1, 0,0, HI- 

Thus, the lattice model provides a self-consistent re- 
alization of the basic assumption underlying engineering 
fracture mechanics, the ability to separate the linear elas- 
ticity calculation from the microscopic physics controlling 
the tip. In general, however, it is hard to test the LEFM 
since it is difficult to reliably measure the SIF in a lattice 
calculation, as this requires an extremely (and imprac- 



tically) large system. Only then is there an appreciable 
range of scales in which one is sufficiently far from the 
crack tip that lattice effects are unimportant and suffi- 
ciently close that the fairly weak square-root singularity 
dominates. However, for the case of a crack accelerating 
from rest, the approach of Kostrov |L| and Eshelby 
provides an analytic prediction of the SIF, independent 
of the details of the dynamics. With this first-principles 
determination of the SIF, it is possible to test the LEFM 
picture. Specifically, we will show that in the lattice 
model, in accord with expectations, the stress intensity 
factor governing the strength of the stress singularity is 
the only information passed from the macroscopic field 
to the process-zone dynamics of the crack tip. In fact, 
this separation is quantitatively accurate even for rather 
small lattices where one might have questioned the effi- 
cacy of the continuum approach. Moreover, these results 
allow us to construct a test of the theory based on the ac- 
tual fracture data presented by Sharon and Fineberg p^] 
for the same case of a crack accelerating from rest. Here, 
however, the data do not appear to conform to theoret- 
ical expectations. At the end, we discuss possible impli- 
cations of this failure. 

We begin with the lattice model. We work with a 2d 
square lattice with unit spacing between the mass points. 
These masses are coupled via both nearest neighbor and 
next-nearest neighbor ideally brittle central force springs 
with spring constants k\ and ki respectively. It is easy to 
show that with the choice k\ = 2/c2 = 2/i = 2A2d, the con- 
tinuum linear elastic limit of this model is isotropic with 
the aforementioned Lame' constants; hence the Poisson 
ratio v^d — 1/4. We will assume that the actual 3d sys- 
tem exhibits plane stress and so can be approximated via 
a 2d system with 
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the actual Poisson ratio of the material being modeled 
is v = V2d/{^ — V2d) = 1/3. The dynamics arising from 
this force is taken to include the possibility of a Kelvin 
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viscosity term. The final equation of motion is therefore 



d 2 u(x) 
dt 2 
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dt 
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((u(x + n) — u{x)) ■ h)h 



((u(x + h) — u{x)) ■ h)n 



(1) 



As discussed in || [lq ], the damping due to nonzero 
Kelvin viscosity occurs only inside the process zone if 
r] is chosen 0(1). Finally, any bond whose length goes 
above a breaking threshold, which we take to be 1, has 
its spring constant set to zero. 

As discussed above, we study in detail a finite length 
crack accelerating from rest. Initially, a crack is placed 
along the mid-line of a sample (here between rows and 
1 of our lattice), extending a length £q from the left edge. 
The top and bottom rows of the lattice have fixed (and 
opposite) displacements u — Ay and the lateral edges are 
free. The loading is chosen to be just below the critical 
value at which the crack will start to propagate; since 
this loading is a decreasing function of crack length, a 
crack with one additional broken bond will in fact start 
to elongate. The system is then allowed to fully relax 
to its equilibrium stress state; this is accomplished via 
a multigrid technique described in detail elsewhere Jl7[. 
Once this is done, an additional bond is broken by hand 
(i.e., its spring constant is set to zero; no actual displace- 
ment of particles is involved) and the crack tip acceler- 
ates. As it moves, we monitor the bonds across the crack 
surface y = 1/2. When all three bonds attached to a 
given point x — (x, 1) break, the crack length is deemed 
to have increased by one and the velocity at that time 
is measured as the inverse of the time interval since the 
last such event. We do not allow bonds off the crack 
line to break, thereby suppressing any possible branch- 
ing instabilities. We also enforce symmetry across the 
crack surface, simulating only the upper half of the lat- 
tice. Each run is characterized by the transverse lattice 
size W, the initial length of the crack £q, damping con- 
stant rj and the driving displacement A. Typical data 
generated by this procedure for both the undamped and 
highly damped cases are presented in Fig. 1. 

According to the classic calculation by Eshelby p[ , 
the stress intensity factor, K, or equivalently the Eshelby 
B factor (which is a constant multiple of K) , at the pro- 
gressing crack tip for a system that starts in equilibrium 
can, up to a certain time, be written as a product of two 
factors: 



Bj(t) = A(v(t))) 
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dx- 



\Jx - £ 



(2) 



Here af™ is the normal stress on the mid-line y = 1/2 
as found from the equilibrium stress field; it diverges of 
course near the edge of the equilibrium crack x = £q 
with a static stress intensity factor. The first factor A(v) 
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FIG. 1: Simulation results for velocity vs. crack length for 
various initial seed cracks. The system widths, W, in lattice 
units are 239, 219, 415, and 1004 for the four runs for de- 
creasing £0. For later convenience, all lengths in the graph 
are scaled by a factor 440 /W. Velocities are normalized to 
the Rayleigh velocity. 



depends only on the instantaneous velocity at time t and 
the second factor (which we will refer to as Bo) depends 
on time only through the instantaneous crack length £(t). 
This equation holds true as long as sound waves reflected 
from the boundaries have not interacted with the crack 
tip. There are two boundaries one must be concerned 
with, the lateral boundaries at y = ±W/2 and the edge 
at x = 0. We take care that all our data come from times 
before these interactions occur. 



Now, if 



the true microscopic breaking events only cou- 



pie to the macroscopic field via Bj, we expect that there 
will be some fixed relationship between the tip velocity 
and this number. Usually, this relationship is considered 
to arise due to energy conservation and is therefore writ- 
ten as an equality of the energy flux into the tip and the 
energy necessary to break bonds along a unit length of 
crack 

I» = f(v)Bj (3) 

Here / is a complicated but explicit function and T is 
the breaking energy. In fact, T can be determined via 
the dependence of steady-state cracks on the driving load 
and thereafter used for the accelerating crack case. For 
our purposes, this hypothesis leads to the existence of a 
universal relationship, independent of W 7 A and £q (but 
dependent on rj) , between the measured velocity v(t) 
and the Eshelby function Bo(£(t)). 

To test this strong prediction of LEFM, we use our 
lattice model simulations discussed above. We calculate 
the integral in Eq. (2) by replacing the continuum stress 
with its lattice analog (defining u = (u, v)) 

a eq( x ) = -hVeq(x,0) ~ [u e q(x - 1, 0)- 

-u eq (x + 1,0)+ v eq {x -1,0) + v eq {x + 1, 0) \4) 

The fields that enter are the equilibrium fields present 
before the crack tip begins to move. The integral is then 
replaced by a sum, taking care to resolve the singulari- 
ties in the integrand. This yields a function Bq(£; W, £q) 
which can be plotted versus the velocity. The results of 
this exercise are presented in Fig. 2, demonstrating al- 
most perfect data collapse even for fairly small transverse 
sizes W. 

So, the lattice model does indeed serve as an instan- 
tiation of the LEFM picture developed by the fracture 
community for straight accelerating cracks. But now, we 
can use the Eshelby Bq function calculated as just de- 
scribed to process actual crack velocity/length data from 
the experiments on PMMA carried out by Sharon and 
Fincbcrg. This data comes from a protocol very similar 
to what we used above in which the loading is set to a 
point where a touch of razor blade will cause the rapid 
elongation of a previously created notch in the material. 
For PMMA, the elastic constants are A = 2680MPa, 
/j = MPa, giving v = .35. Going to the plane stress case 
leads to an effective Poisson ratio V2d — -259 which is 
very close to the value forced upon us by the restriction 
to n.n. and n.n.n. central forces. We did check that gen- 
eralizing our model to include bond-bending springs and 
thereby obtaining the precise value of the Poisson ratio 
did not alter in any way the results which we will now 
present. 

In fact, we designed our simulation to mimic as closely 
as possible the experimental situation. Thus, not only 
were the elastic constants chosen appropriately, the as- 
pect ratio was also fixed to that of the experiment ( 
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FIG. 2: Eshelby data collapse for t] = 0, 1 for the runs shown 
in Fig. 1. Notice the slight systematic dependence on £o for 
the rj — 1 data, presumably arising from the larger process 
zone in this case [HI. 



a width of 440 mm. and a length of 380 mm.). The 
lengths of the initial seed cracks were also chosen to be 
in the same ratio to the width as in the experiment. This 
constraint fixed our choice of the actual width in lattice 
units, as we wanted the initial crack to be large (of order 
10) on the lattice scale. 

Fig. 3 shows our attempt to verify the continuum hy- 
pothesis for the PMMA data. In our opinion, the data 
are quite convincingly inconsistent with the universality 
of the v — Bo relationship. Note that this is the oppo- 
site conclusion from that reached by the experimenters 
themselves, as discussed in Ref. |L5|. There, the exact 
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FIG. 3: Attempted data collapse for the Sharon-Fineberg 
PMMA experiment. Initial crack lengths are given in mm. 



Eshelby function was replaced by an approximate form 
based on the static stress intensity factor for a crack of 
length £(t). Note that these functions while close for 
£ ~ £q, become quite different for larger length cracks. 
The static stress intensity factor saturates for £ ^ W/2, 
whereas the exact Eshelby function increases linearly in 
this regime. We have checked that the use of this approx- 
imation happens to push the data into closer agreement 
with the theory and could lead to a mistaken impres- 
sion of data collapse. This does not happen with the ex- 
act expression and conversely the approximate form does 
not lead to a very good data collapse for the simulation 
study discussed above. The same situation obtains for 
the Sharon-Fineberg measurements on fracture in glass 
(results not shown). There the data collapse is, if any- 
thing, even worse. 

It should be noted that the agreement between the 
lattice theory and LEFM is very dependent on our ab- 
solute suppression of off-axis cracking. If all bonds are 
allowed to break, above a critical velocity the crack no 
longer propagates in a straight line, and the assumptions 
implicit in the Eshelby calculation of the instantaneous 
SIF break down completely. At that point, there is no 
longer any way to reliably measure the SIF and thus test 
LEFM. Needless to say, LEFM in its simplest formu- 
lation does not predict the direction of branching, and 
thus cannot address to the post-instability dynamics of 
the crack. The effects of the instability in the experiment 
(which appear to exhibit different dynamics than in the 
lattice model in this regime) are clearly visible at later 
times (larger velocities), but cannot explain the apparent 
failure of the data collapse at smaller velocities. 



The question remains how to explain this failure of 
LEFM in the Sharon-Fineberg experiment. Putting aside 
for the moment the obvious possibility of systematic er- 
rors in the determination of the velocity or position of the 
crack, the only explanation that suggests itself is that the 
process zone is for some reason very large. The mono- 
tonic dependence on £q, similar to, though much larger 
than that seen in the r/ — 1 simulation, is consistent 
with this hypothesis. This might be due either to the 
polymeric nature of the microstructure of the polymer 
system, the heterogeneity of the small-scale structure, or 
to the large dissipation in the system. However, given 
the success of LEFM on our relatively small systems, it 
is not easy to accept this as a sufficient explanation. Fur- 
ther experiments on other systems are clearly necessary 
for an unraveling of this conundrum. 



The authors wish to thank J. Fineberg and E. Sharon 
for providing the raw data from their experiment, and for 
extensive discussions. The work of DAK is supported in 
part by the Israel Science Foundation. The work of HL is 
supported in part by the NSF, grant no. DMR-0101793. 



[7] 

[8] 

[9] 
[10] 

[11] 



[12] 
[13] 

[14] 
[15] 
[16] 



[17] 



J. Fineberg, S. P. Gross, M. Marder and H. L. Swinney, 
Phys. Rev. Lett. 67, 457 (1991); Phys. Rev. B 45, 5146 
(1992). 

E. Sharon and J. Fineberg, Philos. Mag. 78, 243 (1998). 
J. Fineberg and M. Marder, Phys. Rep. 313, 1 (1999). 
R. D. Deegan, P. J. Petersan, M. Marder, and H. L. Swin- 
ney, Phys. Rev. Lett. 88, 014304 (2002). 
K. B. Broberg, "Cracks and Fracture", (Academic Press, 
1999, San Diego). 

L. I. Slepyan, Doklady Akademii Nauk SSSR 258 561 
(1981) [Sov. Phys. Dokl. 26, 538 (1981)];Sh. A. Ku- 
lamekhtova, V. A. Saraikin, and L. I. Slepyan, Izv. AN 
SSSR. Mekhanika Tverdogo Tela, 19, 112 (1984) [Mech. 
Solids 19, 102 (1984)]. 

M. Marder and S. P. Gross, J. Mech. Phys. Solids 43 1 
(1995). 

D. A. Kessler and H. Levine, Phys. Rev. E 59, 5154 
(1998). 

D. A. Kessler, Phys. Rev. E 61, 2348 (2000). 

L. Pechenik, D. A. Kessler, and H. Levine, Journal of the 
Mechanics and Physics of Solids 50, 583 (2002). 
L. I. Slepyan, M. V. Ayzenburg-Stepanenko, J. P. 
Dempsey, Mechanics of Time-Dependent Materials 3, 
159 (1999). 

E. Gerde and M. Marder, Nature 417, 285 (2001). 

B. V. Kostrov, Prikladnaia Matematika i Mekhanika 30 
1042 (1966) [Applied Math, and Mech. 30, 1241 (1966)]. 
J. D. Eshelby, J. Mech. Phys. Solids 17, 177 (1969). 
E. Sharon and J. Fineberg, Nature 397, 6717 (1999). 
J. S. Langer, Phys. Rev. E 46, 3123 (1992); O. Pla, F. 
Guinea, E. Louis, S. V. Ghasias, and L. M. Sander, Phys. 
Rev. B 57, R13981 (1998). 
D. A. Kessler and H. Levine, in preparation. 



